////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#include "TPCHist.h"
#include "BaseUtils.h"

ClassImp(TPCHist);

Double_t TPCHist::fX0 = -999;
Double_t TPCHist::fX1 = -999;
Double_t TPCHist::fY1 = -999;

TPCHist::TPCHist(const TH2D *horiginal)
{
  const TString pwd = gSystem->pwd();

  if(pwd.Contains("Xmin-9")){
    fX0 = -0.9;
  }
  else if(pwd.Contains("Xmin-8")){
    fX0 = -0.8;
  }
  else if(pwd.Contains("Xmin-5")){
    fX0 = -0.5;
  }
  else if(pwd.Contains("Xmin-3")){
    fX0 = -0.3;
  }
  else if(pwd.Contains("Xmin-2")){
    fX0 = -0.2;
  }
  else if(pwd.Contains("Xmin-1")){
    fX0 = -0.1;
  }
  else if(pwd.Contains("Xmin0")){
    fX0 = 0;
  }
  else if(pwd.Contains("Xmin1")){
    fX0 = 0.1;
  }
  else if(pwd.Contains("Xmin2")){
    fX0 = 0.2;
  }
  else{
    printf("\nTPCHist::TPCHist wrong setting of fX0 %s\n", pwd.Data()); exit(1);
  }

  if(pwd.Contains("Xmax20")){
    fX1 = TMath::Log10(20);//1.3
  }
  else if(pwd.Contains("Xmax40")){
    fX1 = TMath::Log10(40);//1.6
  }
  else{
    //printf("\nTPCHist::TPCHist using default fX1 %e %s\n", fX1, pwd.Data());
    printf("\nTPCHist::TPCHist wrong setting of fX1 %s\n", pwd.Data()); exit(1);
  }

  fY1 = horiginal->GetYaxis()->GetBinUpEdge(horiginal->GetNbinsY())*0.99;

  Double_t anchorN=1e9;
  
  if(pwd.Contains("anchor0")){
    anchorN = 5E2;
  }
  else if(pwd.Contains("anchor1")){
    anchorN = 1E3;
  }
  else if(pwd.Contains("anchor2")) {
    anchorN = 2E3;
  }
  else if(pwd.Contains("anchor3")) {
    anchorN = 5E3;
  }
  else if(pwd.Contains("anchor4")){
    anchorN = 1E4;
  }
  else if(pwd.Contains("anchor5")){
    anchorN = 2E4;
  }
  else{
    printf("\nTPCHist::TPCHist uses default anchorN %e\n", anchorN);
  }

  Int_t ancopt = -999;
  if(pwd.Contains("ancopt0")){
    ancopt = 0;
  }
  else if(pwd.Contains("ancopt1")){
    ancopt = 1;
  }
  else{
    printf("\nTPCHist::TPCHist wrong ancopt!! %s\n", pwd.Data()); exit(1);
  }

  Int_t yrebin = -999;
  if(pwd.Contains("YB1"))
    yrebin = 1;
  else if(pwd.Contains("YB2"))
    yrebin = 2;
  else if(pwd.Contains("YB3"))
    yrebin = 3;
  else if(pwd.Contains("YB4"))
    yrebin = 4;
  else if(pwd.Contains("YB5"))
    yrebin = 5;
  else if(pwd.Contains("YB6"))
    yrebin = 6;
  else{
    printf("\nTPCHist::TPCHist no YB set! %s\n\n", pwd.Data()); exit(1);
  }

  Int_t xrebin = -999;
  if(pwd.Contains("XB1"))
    xrebin = 1;
  else if(pwd.Contains("XB3"))
    xrebin = 3;
  else if(pwd.Contains("XB5"))
    xrebin = 5;
  else if(pwd.Contains("XB7"))
    xrebin = 7;
  else if(pwd.Contains("XB9"))
    xrebin = 9;
  else{
    printf("\nTPCHist::TPCHist no XB set! %s\n\n", pwd.Data()); exit(1);
  }

  Ini(horiginal, xrebin, yrebin, anchorN, ancopt);

  Print();
}

void TPCHist::Ini(const TH2D *horiginal, const Int_t xrebin, const Int_t yrebin, const Double_t anchorN, const Int_t ancopt)
{
 //
  //no modification to horiginal
  //yrebin is to rebin the hist, 
  //{x,y}{0,1} is the range of horiginal to be fitted 
  //anchorN is the weight to regenerate
  //
  //
  
  //highdensity has to be odd so that selected points are equal-distant
  if(xrebin%2==0){
    printf("TPCHist::Ini even highdensity!! The selected points are NOT equal-distant between high- and low-density region! %d\n", xrebin); exit(1);
  }
  if(xrebin<1){
    printf("TPCHist::TPCHist xrebin <=1 %d\n", xrebin);exit(1);
  }

  //not necessary to be divisor, as long as other places is rebins the same.
  //const Int_t nxb = horiginal->GetNbinsX();
  //if(nxb%highdensity!=0){
  //printf("TPCHist::Ini x-rebin not a divisor %d %d\n", nxb, highdensity); exit(1);
  //}

  fHist = horiginal;

  fXrebin = xrebin;
  fYrebin = yrebin;
  fAnchorN = anchorN;
  fAncOpt = ancopt;
  //fLPmax = lpmax;
}

void TPCHist::Print()
{
  printf("\n\nTPCHist::Print()\n"); 
  printf("\n");
  printf("PbinThres()= %e ; \n", PbinThres());
  printf("Xmid()= %e ; \n", Xmid());
  printf("X0()= %e ; \n", X0());
  printf("X1()= %e ; \n", X1());
  printf("Y0()= %e ; \n", Y0());
  printf("Y1()= %e ; \n", Y1());
  printf("\n");
  printf("fXrebin= %d ; \n", fXrebin);
  printf("fYrebin= %d ; \n", fYrebin);
  printf("fAnchorN= %e ; \n", fAnchorN);
  printf("fAncOpt= %d ; \n", fAncOpt);
  printf("\n");
  printf("\n\n");
}


TH2D* TPCHist::GetFitData(TH2D *& hfitdata, TH1D *& hrawyield, TH2D *& hprefilter, TH1D *& hfityield, const Bool_t krest)
{
  //
  //hfitdata is used finally for the fitting, if krest = 1, the entries outside the range are kept in the returned hist
  //hrawyield is the yield before filtering and regenerating
  //hprefilter is the distribution before filtering and regenerating
  //hfityield is the distribution after filtering and regenerating, i.e. the nominal statistics
  //

  hfitdata=(TH2D*)fHist->Clone("hfitdata");

  //==========================================================================
  //rebin raw data
  //==========================================================================
  if(fYrebin>1){
    hfitdata->RebinY(fYrebin);
  }

  /*
  if(fXrebin>1){
    hfitdata->RebinX(fXrebin);
  }

  hrawyield = BaseUtils::GetHyield(hfitdata, "hrawyield");
  hfityield = BaseUtils::Regenerate(hfitdata, fAnchorN,  fAncOpt, PbinThres());

  //dummy
  hprefilter = (TH2D*)fHist->Clone("hprefilter"); hprefilter->Scale(0); 
  TH2D *hrest = (TH2D*)fHist->Clone("hrest"); hrest->Scale(0);

  //test!!! <--
  */

  //_______________________________________________________________
  //fine binnning at low p very important, otherwise pon and electron have much bigger systematic errors  
  //_______________________________________________________________

  Int_t *fitbins=new Int_t[hfitdata->GetNbinsX()];
  const Int_t fitnbs = BaseUtils::RebinX(hfitdata, fXrebin, Xmid(), fitbins);
  if(fitnbs<10){
    printf("TPCHist::GetFitData too few fitnbs %d<10\n", fitnbs); exit(1);
  }

  if(fitnbs>hfitdata->GetNbinsX()){
    printf("TPCHist::GetFitData fitnbs>hfitdata->GetNbinsX() %d %d\n", fitnbs, hfitdata->GetNbinsX()); exit(1);
  }

  for(Int_t ii=0; ii<fitnbs; ii++){
    printf("TPCHist::GetFitData fitbins[%d]= %d ; BinLowEdge: %e ; BinCenter: %e\n", ii, fitbins[ii], hfitdata->GetXaxis()->GetBinLowEdge(fitbins[ii]), hfitdata->GetXaxis()->GetBinCenter(fitbins[ii]));
  }

  //==========================================================================
  //cut out range
  //==========================================================================
  TH2D *hrest = BaseUtils::CutHist(hfitdata, X0(), X1(), Y0(), Y1(), krest);

  //==========================================================================
  //record raw yield
  //==========================================================================
  hrawyield = BaseUtils::GetHyield(hfitdata, "hrawyield");

  //==========================================================================
  //filter out x bins
  //==========================================================================
  hprefilter = (TH2D*)hfitdata->Clone("hfitprefilter");

  BaseUtils::SelectXbins(hfitdata, fitbins);
  delete[] fitbins;
  //==========================================================================
  //regenerate statistics according to fAnchorN
  //==========================================================================
  hfityield = BaseUtils::Regenerate(hfitdata, fAnchorN, fAncOpt, PbinThres());

  //==========================================================================
  //normalize per slice
  //==========================================================================

  //has to normalize after regeneration! because the normalization is not identical to fAnchorN

  TH2D * htmpnor = BaseUtils::NormalHist(hfitdata, PbinThres());
  const TString hfdn = hfitdata->GetName();
  delete hfitdata;
  hfitdata = htmpnor;
  hfitdata->SetName(hfdn);

  return hrest;
}

Bool_t TPCHist::FilterIX(const TH2D * hh, const Int_t ix)
{
  //return 1 if remove, 0 if keep

  //just reject empty bins
  const Double_t cont = hh->Integral(ix, ix, 1, hh->GetNbinsY());

  if(cont < EPSILON)
    return kTRUE;
  else
    return kFALSE;
}
